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Abstract 

This article proposes an adaptive multisensor fusion algorithm for acoustic emission source location in isotropic plate- 
like structures in noisy environments. Overall, the approach consists of the following four main stages: (a) feature extrac- 
tion, (b) sensor selection based on a binary hypothesis testing, (c) sensor weighting based on a well-defined reliability 
function, and (d) estimation of the acoustic emission source location based on the extended Kalman filter. The perfor- 
mance of the proposed algorithm is validated through pencil lead breaks performed on an aluminum plate instrumented 
with a sparse array of piezoelectric sensors. Two experimental setups have been used to simulate a highly noisy environ- 
ment. In the first setup, the experimental signals have been artificially corrupted with different levels of noise generated 
by a Monte Carlo simulation. In the second setup, two piezoelectric transducers have been used to continuously gener- 
ate high-power white noise during testing. The results show the capability of the proposed algorithm for inflight struc- 
tural health monitoring of metallic aircraft panels in highly noisy operational situation. 
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Introduction 

In the past two decades, significant efforts have been 
made toward the development of integrated structural 
health monitoring (SHM) systems in order to increase 
the safety of aircraft systems 1 and reduce life -cycle 
costs. 4 The techniques based on sparse arrays of piezo- 
electric sensors, which have the capability of transmit- 
ting and receiving Lamb waves are among the most 
promising candidates. Damage location using Lamb 
waves can be achieved either by an "active-passive" 
approach or by a "passive-only" approach. In the 
"active-passive" approach, diffractions of piezoelectri- 
cally actuated waves can be used to locate existing dam- 
age in a postmortem mode. 5,6 In the "passive-only" 
approach, growing damage (e.g. fatigue cracks) or sud- 
den impacts can be located by monitoring acoustic 
emissions (AEs) in a real-time mode. 7 15 This article 
focuses on a "passive-only" approach using a sparse 
array of piezoelectric transducers. Traditionally, the 
damage or impact location is based on time-of-flight 
(TOF) triangulation of wave measurements taken at 
multiple receiving points. This technique works very 



well when the wave velocity (V g ) in the test material 
and the arrival time (t,) of the signal at all three sensor 
locations are known. The damage or impact location is 
identified by drawing three circles of radii (R,), whose 
centers coincide with the three sensor locations. The 
radius R, is obtained by multiplying the time of arrival 
of the signal (/,■) with the wave velocity (V g ). The inter- 
section point of these three circles is the damage loca- 
tion (Figure 1(a)). However, the TOF and wave 
velocity are two uncertain parameters. In general, 
uncertainties can be caused by random and systematic 
errors. The random errors are caused by unknown and 
unpredictable changes in the TOF measurements, 
including instrumentation noise, temperature changes, 
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Figure I. AE source location likelihood area: (a) three sensors with different uncertainties and (b) four sensors with different 

uncertainties. 

AE: acoustic emission. 



and so on. 16-21 The systematic errors are mostly caused 
by the digital signal-processing technique used for ana- 
lyzing the time waveforms. 22 ' 23 As a result, rather than 
the damage being located at a single point at the inter- 
section of the circles, it can be located anywhere in the 
dark overlapped region at the intersection of the rings 
around the sensor locations (Figure 1(a)). The ring's 
width represents the uncertainty in the measured dis- 
tance as a result of TOF and/or wave velocity uncer- 
tainties. It should be noted that, uncertainties can be 
directionally dependent on anisotropic structures. In 
general, the lower the signal-to-noise ratio (SNR), the 
higher the ring's width (i.e. higher uncertainty in TOF 
measurements). To overcome these limitations, model- 
based approaches 24 have been proposed for impact 
identification in simple structures. In these methods, the 
location and type of impact in a numerical model are 
iteratively changed until the predicted responses match 
the measured responses. For the damage/impact loca- 
tion in more complex structures, where accurate model- 
based predictions may not be possible, artificial neural 
networks have been proposed. 2526 The time reversal 
approaches have been used to accurately find the source 
location in complex structures. 27-30 Sohn et al. 28 and 
Park et al. 29 have shown that the time reversal algo- 
rithm is able to efficiently localize the impact in isotro- 
pic and anisotropic plate-like structures. However, 
neural networks and time reversal approaches require 
an extensive number of training observations prior to 
deployment, making these methods quite onerous from 
a computational or data storage point of view. Kundu 
et al. proposed an alternative method based on opti- 
mizing an objective function to find the impact location 
in isotropic and anisotropic plates. However, the 



proposed objective function in that reference had the 
inherent problem of multiple singularities, which was 
overcome in Kundu et al. by modifying the objective 
function. The optimization technique was further 
improved by Hajzargarbashi et al. Kundu et al. pro- 
posed a technique to locate acoustic source in large ani- 
sotropic plates that does not require the knowledge of 
the direction-dependent velocity profile or a dense array 
of sensors. In Gaul 34 and Gaul et al., the nonlinear 
least squares optimization adopting the Gauss-Newton 
method was proposed to determine the location, time 
lag, and velocity of "synthetic" AE signals. Ciampa and 
Meo 13 used Newton's iterative method to calculate the 
coordinates of the impact location and the wave velo- 
city. A genetic algorithm (GA) for the optimization of 
an objective function based on the modified triangula- 
tion methodology was proposed by Coverley and 
Staszewski. 36 Dehghan Niri and Salamone 14 proposed 
a probabilistic approach for source localization in the 
isotropic plate-like structures based on the extended 
Kalman filter (EKF) theory. However, they just consid- 
ered the systematic errors in the TOF measurements, 
without taking into account the random errors caused 
by the unknown and unpredictable changes in the mea- 
surements, including environmental noise effect. Here, 
an alternative framework, based on a multisensor data 
fusion algorithm, is proposed for the real-time AE 
source localization and wave velocity estimation in the 
isotropic plate-like structures, in a highly noisy environ- 
ment. The multisensor data fusion is a fairly young 
field, which has mainly been considered and developed 
in military target tracking and autonomous robotics. 37 
Multisensor data fusion is a process similar to the meth- 
ods humans use to integrate data from multiple sources 
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Figure 2. Schematic architecture of proposed adaptive multisensor fusion algorithm. 



and senses to make decisions. In addition to the statisti- 
cal advantage gained by combining same-source data 
(e.g. obtaining an improved estimate of physical phe- 
nomena via redundant observations), 14 the use of multi- 
ple sensors may increase the accuracy with which a 
quantity (e.g. damage location) can be observed and 
characterized (Figure 1(b)). 

In this process, the accuracy of the data is the most 
crucial aspect. In fact, the fusion of inaccurate data, 
such as "noisy" data, can badly affect the estimation of 
the damage location. Therefore, each sensor's informa- 
tion (i.e. TOF measurements) cannot be used equally. 
In fact, the information extracted from each sensor 
depends on the actual SNR, which is related to several 
factors, including intensity of the damage event, envi- 
ronmental and instrumentation noise, attenuation, and 
so on. In the proposed approach, the sensors with 
higher SNR are given more weight, and those with 
smaller SNR are given less weight. For example, in the 
configuration shown in Figure 1(b), the sensor 4 should 
have more weight than the other sensors in the AE 
source location estimation. A schematic diagram of the 
proposed algorithm is shown in Figure 2. 

In particular, the algorithm begins by a feature 
extraction from each recorded signal. These features 
consist of dominant frequency (f,), TOF,, peak ampli- 
tude of the signal envelope (A max ), and standard devia- 
tion of the actual noise level (a,,). Next, a decision is 
made to use or disregard data from a given sensor, 
based on a binary hypothesis test. Then, the informa- 
tion from the remaining sensors is weighted according 
to a well-defined reliability function. Finally, the 
weighted data are used as input to an EKF algorithm, 
which iteratively estimates (a) the AE source location 
and (b) the wave velocity. 

AE source location 

Assuming an arbitrary Cartesian coordinate system, 
the AE source is at unknown coordinates (x s , y s ) in the 
plane of the plate and the sensors are located at (x,-, y,), 



m sensor3 

//(x 3 ,y 3 ) 





sensor5 



Figure 3. AE source location and sparse array of n 
piezoelectric sensors. 
AE: acoustic emission. 

as shown in Figure 3. It is well known that given the 
wave velocity of a particular Lamb wave mode, at least 
three sensors are necessary to locate an AE source in 
the plate-like structures. If t m is the travel time or TOF 
to reach the first sensor (master sensor), and A?„„ are 
the time difference between master sensor and the z'th 
sensor, the following equations can be obtained 



(x m -x s Y + (y m - y s f - (t m V g V = 



(1) 



fc - x s ) 2 + [y t -y s ) 2 - [{t m + Af„„) V g f = (2) 

where V g is the group velocity of the particular Lamb 
wave mode, which is a function of the product of the 
frequency /and the plate thickness d 



V g = F{fd) 



(3) 



By knowing the properties of the plate, mass 
density p, thickness d, Young's modulus E, and 
Poisson's ratio v, V g can be calculated by the Rayleigh- 
Lamb equations* Combining equations (1) and (2), the 
following equation can be derived 



Af„, 



\]{Xi - x s ) 2 + [y t -y s ) 2 - \J{x m -x s ) 2 + (y m -y s ) 2 



(4) 
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In general, the difference between the TOF of the 
master sensor to different sensors (At,,,,) is used to solve 
this set of nonlinear equations with unknowns X = (x s , 
y s , V g ). It should be noted that the master sensor is 
defined as the first sensor triggered by the AE source 
(i.e. the sensor closest to the AE source). However, as 
mentioned in the previous section, several sources of 
error in the TOF measurements (i.e. dispersion, noise, 
and temperature) can affect the accuracy of this solu- 
tion. In this work, to take into account the uncertainty 
in TOF, the arrival time t\ are treated as mutually inde- 
pendent Gaussian random variables. Also considering 
wave velocity and location uncertainties, the unknowns 
x s , t'.s, and V s are initially considered as mutually inde- 
pendent Gaussian random variables. It should be noted 
that the arrival time t\ and the unknowns x s , y s , and V g 
are not independent since those are related by equation 
(4). Based on this assumption, the probability density 
function of the time difference A/„„ can also be defined 
as a Gaussian random variable 41 with mean and var- 
iance defined as 



Af, 1; 



■ti-t„„ er 



A(„, 



2 2 

:0 ^ + < 



(?) 



The determination of time differences At mi and the 
related uncertainties in highly noisy environment is dis- 
cussed in the "Statistical characterization of time of 
flight measurement" section. In this probabilistic frame- 
work, an EKF approach is used for the optimal estima- 
tion of the state vector X= (x s ,y s , V g ). 

Statistical characterization of TOF 
measurement 

The errors in the TOF measurements of Lamb waves 
are strongly dependent on the SNR, as has been quali- 
tatively discussed in the "Introduction" section. In gen- 
eral, the lower the SNR, the higher the uncertainty in 
the TOF measurement. This section presents a statisti- 
cal characterization of TOF errors for different noise 
levels. Toward this end, flexural Lamb waves are first 
modeled analytically to predict the acoustic field in an 
isotropic plate-like structure generated by the AE 
sources. Then, the simulated signals are corrupted arti- 
ficially with different noise levels; using a Monte Carlo 
simulation, the statistical distribution of the TOF error 
for each level of noise is calculated. 

Lamb wave model due to a point source excitation 

In the plate-like structures, the following two funda- 
mental Lamb wave modes can propagate: symmetric 
(S n ) and antisymmetric (A n ). The displacements of the 
symmetric modes occur in the direction of wave propa- 
gation (i.e. extensional modes), while the antisymmetric 



modes have displacements transverse to the wave pro- 
pagation direction (i.e. flexural modes). In this article, 
it is assumed that only the lowest order flexural mode 
(A ) is the dominant mode. 42 The Lamb waves have 
been modeled by Ditri's technique, with normal point 
excitation on the surface of the plate. 43 In particular, to 
model the total out-of-plane surface displacements u- 
associated with an input of finite duration /(/), a 
Fourier transform is used to decompose the input sig- 
nal into its different frequency components, f(co). The 
input, for pencil lead break simulations, has been 
assumed as a one-cycle tone burst enclosed in a 
Hanning window and a center frequency of 1 MHz. 44 
The out-of-plane displacement as a function of time 
and distance from the source is given by 3 



u z [r, t) ■■ 



f[u>) 2^E,„{u>)Hq (k m r) exp {—icot)du> (6) 



where t is the time, r is the distance from the point 
source, k,„ is the wave number of the excited Lamb 
mode (m) corresponding to w, which can be calculated 
through the phase velocity V"' h ( k„, = w/ V"l j , and i is the 
imaginary unit; and E m (w) is the excitability function 
for circular crested waves that depend on the angular 
frequency co and excited mode m. The spatial variation 
in a circular crested Lamb wave model is expressed by 
the Hankel function of first kind H\ . Using Viktorov's 
method, 45 the excitability can be computed. The excit- 
ability of a mode is totally governed by the amount of 
out-of-plane displacement of the mode due to a normal 
point excitation. In this work, in modeling Lamb waves 
generated by pencil lead breaks, it is assumed the A 
mode is the dominant mode. 46 The phase velocity (V p /,) 
for different modes and excitability E(co) just for the 
fundamental mode A , are shown in Figure 4. 

Equation (6) enables the surface displacement to be 
plotted as a function of both distance and time, when 
the plate is excited by a finite duration input at a point. 
Figure 5 shows some of the modeled signals at different 
distances from the source (i.e. pencil lead break). The 
duration and the sampling frequency are fixed to 0.4 ms 
and 3 MHz, respectively. All the signals are normalized 
to the maximum response at 5 cm. These signals will be 
corrupted artificially with different noise levels (gener- 
ated by a Monte Carlo simulation), to determine the 
statistical distribution of the extracted features, includ- 
ing TOF and peak amplitude of the signal envelope. 

TOF measurement 

Several methods have been proposed in the literature 
for the TOF measurement of acoustic signals, including 
threshold crossing, coloration methods, wavelet trans- 
form, and curve fitting approaches. 47 In this work, the 
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mode A in aluminum plate. 
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TOF has been estimated for both experimental and 
numerical signals by the peak amplitude of the signal 
envelope. 36 The envelope can be calculated using the 
definition of analytic signal x A (t) of a real signal x(t) as 
follows 



xa if) =x{t) + ixnit] 



(7) 



where i is the imaginary unit and x H (t) is the Hilbert 
transformation of x(t) defined as 



x H (t) = 



x(t)- 



-dr 



(8) 



A(t). 



x(tf+x H {t) 2 



(9) 



The magnitude of the analytical signal, which is iden- 
tical to the magnitude of the real signal, is called an 
envelope 



The time instant in which the envelope of the signal 
reaches its maximum (A max ), is considered as the TOF. 
Some envelopes of the simulated signals and the corre- 
sponding nominal TOFs are shown in Figure 6. 



TOF error of noisy data 

Several methods have been investigated to analyze the 
error in the TOF measurements of acoustic signals in 
random noise, including threshold crossing 23 and cross 
correlation. 47 ' 48 In this work, the noise is modeled as an 
additive zero-mean Gaussian random process with var- 
iance o- 2 . The ratio between the peak amplitude of the 
signal envelope (A max ) and the standard deviation of 
the noise (a,,) is used as a representative feature of the 
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SNR. In general, the level of the background noise in 
real applications can be estimated by the root mean 
square (RMS) value of a predetermined portion of the 
signal before the first arrival. In this work, the TOF 
error is defined as the difference between the nominal 
TOF captured from predicted signals using equation (6) 
(i.e. noise-free signals) and the TOF extracted from the 
predicted signals artificially corrupted with noise, as 
shown in Figure 7. In particular, a Monte Carlo simula- 
tion has been used to corrupt the predicted signal with 
noise with a standard deviation of 0.01 . 

To capture the statistic distribution of both features 
(i.e. TOF error and (A max /a n )), a predicted signal at a 
given distance r has been corrupted 3000 times with a 
fixed noise level using a Monte Carlo simulation. For 



each simulation, the two features have been extracted. 
Then, histograms have been constructed to predict the 
distribution of the TOF error and A max /cr n . This proce- 
dure has been repeated for different distances (r) and 
different noise levels (o-„). For the sake of clarity, 
Figure 8 shows the envelope of the predicted signal at 
different distances with a superimposed noise level of 
a,, = 0.02; in addition, histograms for some signals 
(i.e. 5, 15, 30, and 50 cm) are indicated. 

It can be observed that as the SNR decreases, the 
distribution of A max ja n for signal plus noise gets closer 
to that of noise alone; as a result, the detection of the 
envelope becomes more unlikely. In addition, as the 
SNR decreases, the standard deviation of the TOF 
error (cr t ) increases. In this work, the standard deviation 
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Figure 9. TOF uncertainty versus distance and noise level. 
TOF: time of flight. 



of the TOF error (o-,) is used as a measure of the TOF 
uncertainty. The variation of cr, versus distance and 
noise level is shown in Figure 9. Note that the further a 
particular sensor is from the AE source, the more unre- 
liable data (TOF) are extracted from it. In the "Sensor 
selection based on binary hypothesis testing" section, a 
binary decision test is proposed to neglect the most 
unreliable sensors. 



Sensor selection based on binary 
hypothesis testing 

Although using more data from a sparse array of sen- 
sors leads to more accurate results, the fusion of inaccu- 
rate data can adversely affect the interpretation of the 
physical phenomenon of interest (i.e. damage/impact 
location). Thus, it would be beneficial in making a deci- 
sion as to whether to use or disregard data provided by 
a given sensor. Toward this end, in this section, a bin- 
ary hypothesis testing is used to classify between a sig- 
nal and noise based on the received signal. The binary 
hypothesis testing are procedures useful for developing 
decision-making algorithms. In particular, these proce- 
dures are used to determine whether a measurement is 
consistent with one of the two states (H and Hi) called 
hypotheses. The hypotheses H and Hi are sometimes 
referred to as the null and the alternative hypotheses, 
respectively. 49 In this context, if v,- is the signal received 
by a given sensor i, the measurements under the two 
different hypotheses are given by 



S + n t 



Hi 
Ho 



(10) 
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where H l denotes the "signal-present" hypothesis and 
H is the null hypothesis (i.e. no signal present); S is an 
arbitrary deterministic signal; and ri\ is zero-mean, 
Gaussian noise with variance o^ . It is assumed that the 
feature (A max /cr n ) extracted from each signal (y,) is a 
sufficient statistic and contains all the required infor- 
mation needed to make a decision. Therefore, the task 
is to decide in favor of H or Hj on the basis of the 
measurements (A max /cr n ). In testing H versus Hi, there 
are two types of errors that can be made as follows: H 
can be falsely rejected or Hj can be falsely rejected. The 
first of these two error types is called a. false alarm (i.e. 
the signal is said to be present when it is not). The sec- 
ond type is called a miss (i.e. the signal is said to be 
absent when it is present). The terms "false alarm" and 
"miss" come from radar problems in which H and Hj 
usually represent the absence and presence of a target, 
respectively. The probability of the first type of error is 
known as the false alarm probability and is denoted as 
P F ; similarly, the probability of the second type of error 
is called the miss probability P M - However, in discuss- 
ing the latter quantity, it is more common to refer to 
the probability of detection, P D = 1 — P M . These 
probabilities can be obtained by the conditional densi- 
ties, P(A,„ ax /a„\ll ) and P{A max ja n \ii\), that reveal how 
the extracted feature (A max /cr„) is distributed under the 
two respective hypotheses. In this work, a Monte Carlo 
simulation has been used to calculate these conditional 
probabilities, considering as deterministic signal (S) the 
Lamb wave modeled as described in the "Lamb wave 
model due to a point source excitation" section for r = 
80 cm. In Figure 10(a), the probability density for noise 
alone (i.e. / > (^„, ua /o"h|Ho)) is plotted along with that for 
signal and noise (i.e. P(A max /o-„\H\)). A threshold A is 
shown. The gray area under the curve for signal plus 



noise represents the P D , while the double-crosshatched 
black area under the curve for noise alone represents 
the P F . Obviously, if the threshold is increased to 
reduce the probability of false alarm, the probability of 
detection will also be reduced. 

The design of a test for H versus Hj involves a 
trade-off between the probabilities of the two types of 
errors since one can be made arbitrarily small at the 
expense of the other. The Neyman-Pearson test criter- 
ion for making this trade-off is to place a bound on the 
false alarm probability, that is 



Pf = 



\Hq \dr = a 



(11) 



Solving this equation for a fixed value of a, the deci- 
sion threshold (A) can be calculated. The decision based 
on the Neyman-Pearson test is to accept H (i.e. 
neglecting the sensor data) if A max /a„ < A or Hj (i.e. 
use the sensor data) if A max /o~„ > A. To describe the 
performance of the decision test, the receiver operating 
characteristics (ROCs), which is nothing more than a 
plot of P D versus P F as threshold A is varied, is shown 
in Figure 10(b) for different levels of noise. It can be 
observed that for a given value of a (i.e. for a given 
threshold A), the level of noise strongly affects the 
detection performance of a given sensor; in particular, 
the higher the noise, the lower the detection perfor- 
mance. This threshold will be fixed in this application 
in the "Results and discussion" section. In addition, it 
should be pointed out that if the number of remaining 
sensors is less than 4, then to solve the set of nonlinear 
equations described in the "AE source location" sec- 
tion (equation (4)), the four sensors that have the 
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largest ratio (A max \(j^) will have to be used. A detailed 
review of the signal detection theory is beyond the 
scope of this article and the readers may refer to 
Trees 49 for more information. 



Sensor weighting 

As mentioned in the previous section, in data fusion 
processes, the accuracy of the data is the most crucial 
aspect. In fact, the fusion of inaccurate data, such as 
"noisy" data, can badly affect the estimation of the 
damage location. Therefore, each sensor's information 
(e.g. TOF measurements) cannot be used equally. In 
this work, it is proposed to achieve an optimal data 
fusion by assigning weighting factors to each sensor's 
information, according to the standard deviation of the 
TOF error (a,). In particular, the sensors with higher 
SNRs (i.e. with smaller a,) are given more weight and 
those with smaller SNRs (with larger a,) are given less 
weight. As a result, the determination of weighting fac- 
tors is related to the estimation of the standard devia- 
tion of the TOF error (a t ) of each sensor. In Figure 8, 
it can be observed that a relationship exists between the 
distribution of the extracted feature (A max /a n ) for signal 
plus noise and a t ; in particular, as the distribution of 
A max /a n for signal plus noise gets closer to that of noise 
alone, the a, becomes larger. Therefore, the mean value 
of the distribution of A,„ ax /a„ can be used for signal 
plus noise, calculated for different noise levels and dif- 
ferent distances, and correlated it with a,. The results 
of this procedure are shown in Figure 1 1 . The exponen- 
tial curve used to fit the data is used to define a reliabil- 
ity function as follows 



a, = 7642e V °» I 



(12) 



Equation (12) is used in the proposed adaptive multi- 
sensor fusion algorithm to weight the contribution of 




Numerical result 
■ Exponential Curve Fitting 



rr,=7642 exp(- 1.069 (A max /<jJ) 



6 7 8 9 10 11 

Mean of A m „Ja„ distribution 



Figure 



Reliability function. 



each data sensor. It should be pointed out that the 
equation (12) is not general, but it depends on the signal 
processing techniques used for the TOF measurements; 
in fact, these techniques have a different sensitivity to 
the noise level. However, in general, a decreasing SNR 
leads to increasing TOF measurement uncertainty; 50 
this fact is clearly considered in this equation. Thus, 
when there is no such a function for the specific signal 
processing and configuration, using this equation as a 
primary function for weighting each sensor data in a 
highly noisy environment is recommended because it 
changes the contribution of each sensor relative to the 
other sensors (its absolute value is not as important as 
its relative value). In particular, the proposed approach 
to perform a multisensor data fusion is based on the 
EKF theory. 



EKF 

In general, a Kalman filter (KF) is a recursive data pro- 
cessing algorithm that estimates the state of a noisy 
dynamic system. 51 In this article, the state of a system 
is used to mean a vector X consisting of n variables 
describing some interesting properties of the system. To 
estimate the state, a KF processes all available mea- 
surements, both accurate and inaccurate. It is note- 
worthy that the KF belongs to the Bayesian inference- 
based estimators that have been widely used in SHM. 52 
The KF has been extensively used in a wide range of 
applications due to its simplicity, optimality, robust- 
ness, sensor fusion ability, and capability of filtering 
out the uncertainties due to the mathematical model 
and measurements. 51,53 55 In general, the KF is a two- 
step process as follows: (a) state prediction according to 
a mathematical model and (b) correction of the state 
according to the measurements collected by the sensors. 
More specifically, in the prediction step, the KF esti- 
mates the state of the system at a given time instant 
and then obtains a feedback control in the correction 
step, by incorporating a new measurement result into 
the a priori estimate to gain an improved a posteriori 
estimate. Although the underlying approach is very 
promising, there is a critical limitation. In fact, the KF 
assumes that the system and measurements are ade- 
quately modeled by a linear dynamic system. The EKF 
is similar to the KF, but it can be used in nonlinear sys- 
tems with linearization using first-order Taylor series 
expansion. 51 For this reason, the EKF has been used to 
handle the set of nonlinear equations (see equation (4)). 
In this work, the state of the system X consists of three 
parameters, namely, the AE source location (x s , y s ) and 
the wave velocity V g . Given a sparse array of N sen- 
sors, the time differences A?,,,,- between the z'th sensor 
and the master sensor represent the measurement data 
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vector Z= [At mi ] {N _ l) , N > 4. The state of the system X 
is related to the measurement vector Z through the 
nonlinear measurement function /; defined in equation (4). 
The key idea underlying the proposed approach is that 
the AE source location can be considered as an unmo- 
vable (static) or slowly fluctuating point. Under this 
assumption, the EKF algorithm can be simplified as 
follows 14 ' 15 



Z k = h(X k ^) 

Xk = %k~ 1 + Kk (Z — Zk 

P k =(l-K k H k )P k _ 1 
dh(X) 



H, 



ax 



(13) 
(14) 
(15) 

(16) 



x=x t - 



Z k is the predicted measurement vector in the Arth 
step according to state estimation X^ \ in the (k — l)th 
step using the nonlinear measurement function (h) 
defined in equation (4); this equation predicts difference 
in the TOFs or components of measurement vector 
based on the estimated state vector including location 
and wave velocity. X k is the estimated states in the /rth 
step derived by correcting the estimated states at the 
(k — l)th step according to the difference between the 
actual measurement Z and the predicted measurement 
vector Z k , and Kalman gain K k is defined as follows 



K k =P k ,H H k P, 



k-\n k 



ikr k - 



,H' k+ R 



(17) 



where P k is the error co variance matrix in the k\h step. 
The Jacobian matrix H k contains the partial derivatives 
of the measurement function h with respect to the 
state X. The objective is to estimate the state vector 
X= [x s , y s , Vg\ given the measurement vector Z and 
covariance matrix R. Note that the measurement vector 
Zin equation ( 14) does not have the subscript k because 
its (N — 1) components (difference in arrival times) can 
be considered constant once the AE signals have been 
acquired. The uncertainty of the measurement vector Z 
is modeled as a zero-mean white Gaussian noise with 
covariance matrix R with diagonal terms defined in 
equation (5) 



R 



'At. 











'At, 



(18) 



As discussed in the "Sensor weighting" section, a ti 
in equation (5) can be calculated using the reliability 
function defined in equation (12). Note that the 
Kalman gain in equation (17) is inversely related to R. 
Therefore, the sensors with a higher SNR (i.e. with 



smaller <r t ) are given more weight, and those with a 
smaller SNR (with larger o-,) are given less weight. The 
EKF will iteratively correct a priori knowledge of the 
state vector X and covariance matrix P, with respect to 
Z and R. The iteration begins with the initialization of 
the state vector estimate X Q and its covariance matrix 
P as follows 



-^o- [x So ,y So , V g0 ] 

a k ° 

0° oi 



■ v Vi 



(19a) 



(19b) 



It should be observed that for large values of 
A max /cr n , the probability of detection approximately 
tends to 1 and a, converge to zero; therefore, singulari- 
ties may occur in the EKF algorithm. To overcome this 
limitation, it is proposed in this work to fix a,, = 1 |xs 
as A max /(T n > 8. In fact, in the case of a sufficiently high 
SNR, each sensor's data can be treated equally. 15 

Initiation of location 

The starting estimates of the AE source location x So and 
y So in equation (19a) can be calculated by geometrical 
considerations. Given an arbitrary group of N sensors 
at locations (Xi,yt), it can be assumed that the estimated 
coordinates of the source location x s and y s can take 
any value in the intervals of (x sL , x sU ) and {y sL , y sU ), 
where x sL and y sL are the lower bounds defined as 

x sL - min ( [xi] N ) ; y sL -min{ \y t ] N ) (20) 

and x sU and y s y are the upper bounds defined as 

x sU -max([xi] N );y sU -max(\yi] N ) (21) 

Under these assumptions, the mean and variance for 
an uncorrelated uniform distribution can be expressed 

as follows 



„ _ (xsu+x s l) „ _ (ysu +y S L) fll , 

X S - , » ys a - - \J-!-) 



2 1 / \2 2 1 / ^2 

% = ^2 ( XsU ~ XsL > ' °>'>-o = n ^ sU ~ ysL > 



(23) 



Initiation of group velocity 

The starting estimates of the group wave velocity V g0 
and its variance o>^ j n equation (19b) can be calcu- 
lated by equation (3). In particular, assuming/ as the 
frequency corresponding to the maximum amplitude of 
the fast fourier transform (FFT) of the signal received 
(y,), its velocity V gj can be obtained from the dispersion 
curves defined by equation (3) as follows 
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Figure 1 2. (a) Time history of recorded AE signals for eight sensors without artificial additive Gaussian noise, (b) signals envelope 
using Hilbert Transform, (c) FFT of recorded signals and the dominant frequencies, and (d) initiation of group wave velocity of Lamb 
mode A . 



V g =F{f i d) 



Therefore, given a sparse array of N sensors, the ini- 
tial group wave velocity V gli is calculated as the average 
of the velocity (V gj ), that is 



V - V - 



N 

_i 

N 



and the variance a 2 ,, as 



N 

i 



f> \2 



v g ) 



N 



(24) For the sake of clarity, Figures 12 shows this proce- 
dure for an array of eight sensors. 

In particular, Figure 12(a) shows the time wave- 
forms generated by a pencil lead break on the surface 
of an aluminum plate. Figure 12(b) shows the envelope 
for each signal, and in Figure 12(c) and (d), the FFT 
and the corresponding velocity (V gl ) through the disper- 
sion curves are depicted. In this work, the fundamental 

(25) antisymmetric mode (A ) has been assumed as the 
dominant Lamb wave mode. Overall, the proposed 
approach consists of four main stages as shown in the 
flowchart in Figure 13. The initial wave velocity distri- 
bution is just an initial guess to begin the algorithm. 
Indeed, the wave velocity of propagating waves will be 

(26) estimated at the end of the iteration. The estimated 
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Figure 1 3. Proposed approach flowchart. 
EKF: extended Kalman filter; TOF: time of flight. 
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wave velocity is the average velocity of propagating 
wave packet to reach each sensor that corresponds to 
the related measured TOF of the sensor. 

The first stage (i.e. feature extraction) provides infor- 
mation about the TOF,, dominant frequencies (f t ), peak 
amplitude of the signal envelope (A,„ ax ), and noise level 
(<r„); during the second and third stages (i.e. sensor 
selection and sensor weighting), these features are used 
to eliminate and/or weight the sensor's data as dis- 
cussed in the "Sensor selection based on binary hypoth- 
esis testing" and "Sensor weighting" sections; finally, 
the last stage of the process (i.e. estimation) uses the 
remaining sensor data and a priori information about 
the states, to iteratively estimate the state vector accord- 
ing to the measured information obtained during the 
second stage. The recursive procedure stops when a 
suitable finish condition is met; the elements of the state 
vector X corresponding to the last state estimate give 
the desired parameters of the AE location. In particu- 
lar, to set the finish criterion, two approaches can be 
used as follows: (a) a minimum number of iterations 
can be fixed or (b) a value can be defined as the modu- 
lus of the difference between the state vector estimates 
provided by two consecutive Kalman filtering loops, 
that is 



Finish Criterion = \\X/ S — Xt-\ (27) 

This value is compared to a proper threshold value 
(empirical tests suggest a value of 10~ 5 ). If the 



difference is lower than the threshold, the recursive 
procedure stops, and the best estimate of the three 
parameters is delivered; otherwise, a new Kalman filter- 
ing loop is executed. In equation (27), | denotes the 
norm one. 



Experimental setup 

The experiments were carried out using an aluminum 
plate with dimensions of 910 mm X 910 mm X 3.175 
mm to validate the proposed algorithm. Two different 
experimental setups were used. In the first setup, an 
array of eight piezoelectric transducers in a square con- 
figuration with 560 mm dimension was used, as shown 
in Figure 14. In the second setup, six sensors in hexago- 
nal configuration with the edge of 300 mm were used, 
as shown in Figure 15. Each sensor was connected to a 
preamplifier. For the data acquisition, an eight-channel 
high-speed data acquisition board (Physical Acoustics 
Corporation Micro-II PAC) with a sampling frequency 
of 3 MHz and dedicated software for signal processing 
and storage (AE win ) were used. During testing, the AE 
sources were generated by pencil lead breaks at sys- 
tematic grid locations. To preferentially generate the 
antisymmetric (A ) Lamb wave mode, the lead was 
fractured on the plate surface. 42 Postprocessing of the 
received signals was performed with a personal com- 
puter (PC) running a MATLAB software code imple- 
mented by the authors. Figures 14(b) and 15(b) show 
the schematic experimental setup. In these figures, the 
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Figure 14. (a) The first experimental setup and (b) actual locations of simulated AE source with pencil lead break. 
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Figure 1 5. (a) The second experimental setup and (b) actual locations of simulated AE source with pencil lead break. 



actual source location is marked with an asterisk (*) 
and each sensor is marked with a circle (o). In the first 
setup (Figure 14), the lower bounds x sL and y sL are the 
location of sensor 1, and the upper bounds x sU and y sU 
are the location of sensor 8. In the second experimental 
setup (Figure 15), the lower bounds x sL and y sL are the 
location of sensors 3 and 1, respectively, and the upper 
bounds x s u an d y.iu are tne location of sensors 4 and 6, 
respectively. 

The noise in the first experimental setup was simu- 
lated by artificially corrupting the acquired signals with 
additive Gaussian random noise. However, to simulate 
more realistic situation, in the second setup, two 
sources of noise were applied using two piezoelectric 



transducers generating continuously high level of white 
noise by a two-channel arbitrary waveform generator 
during the pencil lead breaks, as shown in Figure 15. 
The bandwidths of the piezoelectric transducers for 
noise generation were 60 and 1 50 kHz and were used as 
the first and second noise sources, respectively. 



Results and discussion 

In both experimental setups, the proposed algorithm 
was applied in two modes as follows: (a) regular multi- 
sensor fusion and (b) adaptive multisensor fusion . In the 
first mode (i.e. regular multisensor fusion), each 
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sensor's information has been equally used (i.e. a, was 
fixed to 1 |jls). 15 In the second mode (i.e. adaptive multi- 
sensor fusion), the probability of false alarm (P F ) in the 
selection stage was fixed to a = 0.04, which corre- 
sponds to a threshold A = 4 (i.e. solving equation (11)). 



Artificial noise simulation (first experimental setup) 

To validate the proposed algorithm, three levels of 
noise a„ = 0.03, 0.04, and 0.05 were considered. Each 
generated AE signal (i.e. generated by a pencil lead 
break) was corrupted by an additive zero-mean 
Gaussian noise with one noise level using a Monte 
Carlo simulation. In particular, 1000 simulations were 
conducted for each actual location, and the estimated 
locations using the proposed algorithm were plotted. 
Figure 16 shows the estimated AE source locations 
using the "regular multisensor fusion algorithm" and 
the "adaptive multisensor fusion algorithm." The 
improvement obtained as a result of the adaptive pro- 
cedure for all points and for different noise levels is 
clearly observed. Note that large errors were observed 
in some points (i.e. 1, 4, 13, and 16) due to multiple 
reflections from the plate edges. 

To evaluate this algorithm quantitatively, an error s n 
is defined for each estimated point as follows 



£,., = yj (x s - x Sl f + (y s - y s ) i=\ ,N S (28) 

where N s is the number of simulations. The statistics of 
the error £,., mean and standard deviation, are defined 
as follows 



N, 
E £ r, 



— 1 
£,: 



N s 



N s 



E Or, - £,-) 
^=^ 



(29) 



(30) 



The most desirable estimator is one with the lowest 
mean (unbiased) and standard deviation of error. 
Figure 17 shows the mean and standard deviation 
defined in equations (29) and (30) for the two modes 
and noise levels considered. It can be observed that 
good improvements are achieved in terms of the mean 
and standard deviation of estimation error at any 
point. The computational cost of the proposed algo- 
rithm, using a regular PC, is less than 0.03 s for each 
point. Table 1 summarizes the average of the mean and 
standard deviation of radius of error for 16 points for 
each noise level. 



Real noise simulation (second experimental setup) 

To validate the proposed algorithm in a more realistic 
noisy environment, two piezoelectric transducers were 
used in the second experimental setup (see Figure 15) to 
continuously generate two white noise signals by a two- 
channel arbitrary waveform generator. The amplitude 
of the two signals was 0.1-V peak to peak. The transdu- 
cers were connected to two 40-dB amplifiers. The gener- 
ated noise bandwidth of the first (60 kHz bandwidth 
transducer) and second (150 kHz bandwidth transducer) 
sources was set to 20 and 200 kHz, respectively, to simu- 
late both high and low bandwidth noises. Ten pencil 
lead breaks were conducted (N s = 10) on each point, 
and the signals were stored in a PC for postprocessing. 
It should be mentioned that using an artificial noise, the 
information of the noise level is a priori knowledge for 
each sensor; however, in a real situation, the level of 
noise has to be measured for each sensor. The level of 
the background noise a„ for each sensor was measured 
by calculating the RMS value of a predetermined por- 
tion of the signal before the first arrival of each signal, 
so that obviously the level of the noise for each sensor is 
different. Figure 18 shows the estimated AE source loca- 
tions using the "regular multisensor fusion algorithm" 
and the "adaptive multisensor fusion algorithm." The 
improvement obtained as a result of the adaptive proce- 
dure for all points in real noisy environment is clearly 
observed (note that for clarity in Figure 18, each esti- 
mated location is connected with a line to the actual cor- 
responding source location). Figure 19 shows the mean 
and standard deviation defined in equations (29) and 
(30) in the two modes for real noise simulation. The sig- 
nificant error reduction in terms of mean and standard 
deviation of radius of error was achieved. 

Table 2 summarizes the average of the mean and stan- 
dard deviation of radius of error for nine-point real noise 
simulation. It is worth noting that in the second experi- 
mental setup due to high level of noise, a commercially 
available AE software could not even properly operate 
as a result of many false alarms or miss detections. 



Conclusion 

In this article, an adaptive multisensor fusion algorithm 
for locating the AE sources in noisy operational condi- 
tions is proposed. Based on a fixed probability of false 
alarm, a threshold has been calculated to eliminate the 
data from unreliable sensors using the Neyman-Pearson 
test. Then, the remaining data sensors were weighted 
based on a predetermined reliability function that cor- 
relates the ratio between the peak amplitude of the sig- 
nal envelope and noise level, to uncertainty of the TOF 
measurement. Finally, the remaining weighted data are 
fed into an estimation procedure based on the EKF. 
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Figure 16. Actual and estimated locations for 16 simulated AE sources for 1000 Monte Carlo simulations using: (a) noise level I 
and regular sensor fusion, (b) noise level I and adaptive sensor fusion, (c) noise level 2 and regular sensor fusion, (d) noise level 2 
and adaptive sensor fusion, (e) noise level 3 and regular sensor fusion, and (f) noise level 3 and adaptive sensor fusion. 
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Figure 1 7. Mean and standard deviation of error e r for regular sensor fusion and adaptive sensor fusion: (a) noise level I and mean 
of £ r ; (b) noise level I and standard deviation of s r ; (c) noise level 2 and mean of e r ; (d) noise level 2 and standard deviation of e, r ; 
(e) noise level 3 and mean of s r ; and (f) noise level 3 and standard deviation of s r 



The main advantages of the proposed algorithm are 
given as follows: (a) flexibility in fusion of the data due 
to the matrix structure of the estimation stage, (b) 



accurate AE source localization in noisy environments, 
and (c) the computation cost makes it feasible for real- 
time applications. 
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Figure 18. Actual and estimated locations for nine simulated AE sources for 10 pencil lead breaks and real noise simulations using 
(a) regular sensor fusion and (b) adaptive sensor fusion. 
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